## Replication Code for: "The Financial Consequences of Rating International Institutions: Competition, Collaboration, and the Politics of Assessment"
## Ranjit Lall
## ISQ
## October 2020

## Note: figure PDF files will be saved in the working directory; tables are based on the stargazer() console output

## Load R libraries
library(foreign)
library(sandwich)
library(multiwayvcov)
library(arm)
library(dotwhisker)
library(lmtest)
library(Hmisc)
library(ggplot2)
library(gridExtra)
library(plyr)
library(directlabels)
library(matrixStats)
library(dplyr)
library(broom)
library(rdd)
library(tidyr)
library(RColorBrewer)
library(interflex)
library(stargazer)

## Set working directory
setwd()

## Ratings
ratings <- read.csv("ratings_data.csv", header=T)
head(ratings)
dim(ratings)

## Generate averages
## SWE
ratings$swe_mean <- rowMeans(ratings[,c("swe_int", "swe_ext")], na.rm = T)

## MOS
ratings$mos_mean <- rowMeans(ratings[,c("mos_PDR", "mos_CTFR", "mos_FTP", "mos_RAD",  "mos_FA", "mos_UPI", "mos_AP", "mos_CPD", "mos_HP", "mos_EXR", "mos_PPI", "mos_DLL")], na.rm = T)

## MOR
ratings$mor_mean <- rowMeans(ratings[,c("mor_PDR", "mor_FTP", "mor_CPFR", "mor_RAD", "mor_FA", "mor_UPI", "mor_MHR", "mor_EXR", "mor_PPI", "mor_DLL")], na.rm = T)

## Intercorrelations
sum_ind <- c("uk_vfm", "aus_vfm", "den_eff", "net_eff", "swe_mean", "mos_mean", "mor_mean")
df <- as.matrix(ratings[sum_ind])
cor.mat <- rcorr(df)$r

## Mean
result_r <- vector("numeric", 7)
result_p <- vector("numeric", 7)
for(i in 1:length(sum_ind)){
	cor <- mean(cor.mat[-i,i])
	result_r[i] <- cor
}

## P-values
cor.pval.mat <- rcorr(df)$P
sum(cor.pval.mat>=0.1, na.rm=T)/2

## Mean
for(i in 1:length(sum_ind)){
	pval <- mean(cor.pval.mat[-i,i])
	result_p[i] <- pval
}
rbind(result_r, result_p)

## Within-assessment correlations
mos <- as.matrix(ratings[,c("mos_PDR", "mos_CTFR", "mos_FTP", "mos_RAD",  "mos_FA", "mos_UPI", "mos_AP", "mos_CPD", "mos_HP", "mos_EXR", "mos_PPI", "mos_DLL")])
mor <- as.matrix(ratings[,c("mor_PDR", "mor_FTP", "mor_CPFR", "mor_RAD", "mor_FA", "mor_UPI", "mor_MHR", "mor_EXR", "mor_PPI", "mor_DLL")])

## SWE
mean_swe <- cor(ratings$swe_int, ratings$swe_ext, use = "complete.obs")
p_swe <- cor.test(ratings$swe_int, ratings$swe_ext)$p.val
## 1 out of 1 SS

## MOS
mean_mos <- mean(rcorr(mos)$r)
mos_mat <- rcorr(mos)$P
sum(mos_mat>=0.1, na.rm=T)/2
sum(is.na(mos_mat)==F)/2
## 66/66 are SS

## MOR
mean_mor <- mean(rcorr(mor)$r)
mor_mat <- rcorr(mor)$P
sum(mor_mat>=0.1, na.rm=T)/2
sum(is.na(mor_mat)==F)/2
## 37/45 are SS

## Standardize indicators
sum_ind <- c("uk_vfm", "aus_vfm", "den_eff", "net_eff", "swe_mean", "mos_mean", "mor_mean")
sum_ind <- scale(ratings[sum_ind])
ratings$avg_measure <- rowMeans(sum_ind, na.rm = T)

data <- read.csv("indvar_data.csv", header=T)
data <- merge(data, ratings, all = T, by=c("IN", "year"))
data$ln_cont_ia <- log(data$cont_ia+1)
post <- data %>% group_by(IN) %>% filter(is.na(avg_measure)==FALSE) %>% slice(2:n())
pre_5 <- data%>% group_by(IN) %>% filter(is.na(avg_measure)==TRUE) %>% do(tail(., n=5))
post <- aggregate(post[, "ln_cont_ia"], list(IN= post$IN), mean, na.rm = T)
pre_5 <- aggregate(pre_5[, "ln_cont_ia"], list(IN= pre_5$IN), mean, na.rm = T)
avg_measure <- aggregate(ratings[, "avg_measure"], list(IN= ratings$IN), mean, na.rm = T)
post <- plyr::rename(post, c("ln_cont_ia"="post_avg"))
pre_5 <- plyr::rename(pre_5, c("ln_cont_ia"="pre_5_avg"))
avg_measure <- plyr::rename(avg_measure, c("x"="avg_measure"))

fig1 <- merge(post, pre_5, by="IN")
fig1 <- merge(fig1, avg_measure, by="IN")
fig1$avg_d5 <- (fig1$post_avg/fig1$pre_5_avg)-1

## Subset
subset <- fig1[-c(22, 11),]
cor.test(subset$avg_measure, subset$avg_d5)

## FIGURE 1
pdf("Figure 1.pdf")
ggplot(subset, aes(x= avg_measure, y= avg_d5)) + geom_point(shape=20, size = 0) + labs(x = "Average performance rating (standardized)", y = "Average post-assessment change in log contributions") + theme_classic() + theme(axis.line=element_line(size=0.3)) + geom_dl(aes(label=IN), method = list("top.bumpup", cex = 0.65, vjust = -0)) + geom_hline(yintercept=0, linetype = "dashed", size = 0.3, col = "darkgray") + geom_vline(xintercept=0, linetype = "dashed", size = 0.3, col = "darkgray") + geom_smooth(method=lm, size = 0.2, lty = 1, col = "black", level = 0.95, fill = "lightgray")
dev.off()

## Subset ratings
ins <- unique(data$IN)
for(i in 1:length(ins)){
	trends <- paste("trends_", ins[i], sep="")
	assign(trends, subset(data, IN == ins[i]))
}

## FIGURE 2
pdf("Figure 2.pdf", width=12, height=16)
par(mfrow=c(4,3), oma=c(2.5,0,0,0))
plot(trends_UNHCR$year, trends_UNHCR$cont_ia, pch = 17, xlab = "Year", ylab = "Contributions (inflation-adjusted USD millions)", ylim=c(500, 3000), main = "UNHCR: High average rating", col = "black")
lines(trends_UNHCR$year, trends_UNHCR$cont_ia, col = "black")
abline(v=2008, lty=2, col = "gray50")
plot(trends_EBRD$year, trends_EBRD$cont_ia, pch = 17, xlab = "Year", ylab = NA, ylim=c(0, 700), main = "EBRD: High average rating", col = "black")
lines(trends_EBRD$year, trends_EBRD$cont_ia, col = "black")
abline(v=2008, lty=2, col = "gray50")
plot(trends_WB$year, trends_WB$cont_ia, pch = 17, xlab="Year", ylab=NA, ylim=c(3000, 8000), main = "WB: High average rating", col = "black")
lines(trends_WB$year, trends_WB$cont_ia, col = "black")
abline(v=2008, lty=2, col = "gray50")
plot(trends_UNAIDS$year, trends_UNAIDS$cont_ia, pch = 17, xlab ="Year", ylab = "Contributions (inflation-adjusted USD millions)", ylim=c(50, 300), main = "UNAIDS: Low average rating", col = "black")
lines(trends_UNAIDS$year, trends_UNAIDS$cont_ia, col = "black")
abline(v=2008, lty=2, col = "gray50")
plot(trends_COMSEC$year, trends_COMSEC$cont_ia, pch = 17, xlab="Year", ylab=NA, ylim=c(30, 80), xlim=c(2001, 2016), main = "COMSEC: Low average rating", col = "black")
lines(trends_COMSEC$year, trends_COMSEC$cont_ia, col = "black")
abline(v=2011, lty=2, col = "gray50")
plot(trends_UNDP$year, trends_UNDP$cont_ia, pch = 17, xlab="Year", ylab = NA, ylim=c(2000, 5000), main = "UNDP: Low average rating", col = "black")
lines(trends_UNDP$year, trends_UNDP$cont_ia, col = "black")
abline(v=2008, lty=2, col = "gray50")
plot(trends_EDF$year, trends_EDF$cont_ia, pch = 17, xlab="Year", ylab = "Contributions (inflation-adjusted USD millions)", ylim=c(500, 4500), main = "EDF: High average rating", col = "darkgray")
lines(trends_EDF$year, trends_EDF$cont_ia, col = "darkgray")
abline(v=2011, lty=2, col = "gray50")
plot(trends_MLF$year, trends_MLF$cont_ia, pch = 17, xlab="Year", ylab = NA, xlim=c(2006, 2016), ylim=c(50, 150), main = "MLF: High average rating", col = "darkgray")
lines(trends_MLF$year, trends_MLF$cont_ia, col = "darkgray")
abline(v=2012, lty=2, col = "gray50")
plot(trends_IFRC$year, trends_IFRC$cont_ia, pch = 17, xlab="Year", ylab=NA, ylim=c(0, 600), main = "IFRC: High average rating", col = "darkgray")
lines(trends_IFRC$year, trends_IFRC$cont_ia, col = "darkgray")
abline(v=2011, lty=2, col = "gray50")
plot(trends_FAO$year, trends_FAO$cont_ia, pch = 17, xlab="Year", ylab="Contributions (inflation-adjusted USD millions)", ylim=c(400, 1400), main = "FAO: Low average rating", col = "darkgray")
lines(trends_FAO$year, trends_FAO$cont_ia, col = "darkgray")
abline(v=2008, lty=2, col = "gray50")
plot(trends_UNEP$year, trends_UNEP$cont_ia, pch = 17, xlab="Year", ylab=NA, ylim=c(100, 600), main = "UNEP: Low average rating", col = "darkgray")
lines(trends_UNEP$year, trends_UNEP$cont_ia, col = "darkgray")
abline(v=2008, lty=2, col = "gray50")
plot(trends_UNW$year, trends_UNW$cont_ia, pch = 17, xlab="Year", ylab = NA, ylim=c(0, 300), main = "UN Women: Low average rating", col = "darkgray")
lines(trends_UNW$year, trends_UNW$cont_ia, col = "darkgray")
abline(v=2008, lty=2, col = "gray50")
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend("bottom", legend="First assessment", lty=2, horiz=T, cex=1.5, xpd=T, bty="n", col="gray50")
dev.off()

## Normalize alliances
stdize = function(x, ...) {(x - min(x, ...)) / (max(x, ...) - min(x, ...))}
data$alliances <- stdize(data$alliances,na.rm=T)

## Check correlation
cor.test(data$alliances, data$competition)

## Define samples
assess <- c("uk", "aus", "den", "net", "swe", "mos", "mor", "data")
perf_vars <- c("uk_vfm", "aus_vfm", "den_eff", "net_eff", "swe_mean", "mos_mean", "mor_mean", "avg_measure")

## Create treatments
for(i in 1:length(perf_vars)){
	subset <- data %>% group_by(IN) %>% filter(sum(get(perf_vars[i]), na.rm=T)!=0)
	subset[[paste("rating")]] <- ifelse(is.na(subset[[perf_vars[i]]])==F, scale(subset[[perf_vars[i]]]), 0)
	assign(paste(assess[i]), subset)
}

## Create leads and lags
for(i in 1:length(assess)){
	for(j in 1:10){
	rating_lg <- paste("rating_lg", j, sep="")	
	rating_ld <- paste("rating_ld", j, sep="")
	lls <- get(assess[i]) %>% group_by(IN) %>% mutate(!!rating_lg := lag(rating, order_by=IN, n=j), !!rating_ld := lead(rating, order_by=IN, n=j))
	assign(paste(assess[i]), lls)
		}
}

## Effect of moderators on treatment
## Lag moderators and outcome
for(i in 1:length(assess)){
	comp_lg <- paste("comp_lg", sep="")	
	all_lg <- paste("all_lg", sep="")
	ln_cont_ia_lg <- paste("ln_cont_ia_lg", sep="")
	lags <- get(assess[i]) %>% group_by(IN) %>% mutate(!!comp_lg := lag(competition, order_by=IN), !!all_lg := lag(alliances, order_by=IN), !!ln_cont_ia_lg := lag(ln_cont_ia, order_by=IN))
	assign(paste(assess[i]), lags)
}

## TABLE A5 (Appendix)
mts.lm <- mts.coef <- mts.se <- mts.p <- list()
for(i in 1:length(assess)){
	mts.lm[[i]] <- lm(as.formula(paste("rating~comp_lg+all_lg+factor(IN)+factor(year)")), data=get(assess[i]))
	mts.coef[[i]] <- coeftest(mts.lm[[i]], cluster.vcov(mts.lm[[i]], get(assess[i])$IN))
	mts.se[[i]] <- mts.coef[[i]][,2]
	mts.p[[i]] <- mts.coef[[i]][,4]
}

stargazer(mts.lm, se=mts.se, p=mts.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Effect of outcome on treatment assignment 
## TABLE A6 (Appendix)
ota.lm <- ota.coef <- ota.se <- ota.p <- list()
for(i in 1:length(assess)){
	ota.lm[[i]] <- lm(as.formula(paste("rating~ln_cont_ia_lg+factor(IN)+factor(year)")), data=get(assess[i]))
	ota.coef[[i]] <- coeftest(ota.lm[[i]], cluster.vcov(ota.lm[[i]], get(assess[i])$IN))
	ota.se[[i]] <- ota.coef[[i]][,2]
	ota.p[[i]] <- ota.coef[[i]][,4]
}

stargazer(ota.lm, se = ota.se, p=ota.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Baseline models
## Eq. 3
base3.lm <- base3.coef <- list()
for(i in 1:length(assess)){
	base3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)")), data=get(assess[i]))
	base3.coef[[i]] <- tidy(coeftest(base3.lm[[i]], cluster.vcov(base3.lm[[i]], get(assess[i])$IN)))
}

## Coefficient plot
l <- lapply(base3.coef, function(x) {i <- which(x$term == "rating_lg1")
	if (length(i) > 0) x[i, ]
	else NA })
rating_lgs <- rbind(l[[1]], l[[2]], l[[3]], l[[4]], l[[5]], l[[6]], l[[7]], l[[8]])
rating_lgs$term <- assess

coefplot1 <- dwplot(rating_lgs) + theme(panel.background = element_rect(fill = NA, colour = "black")) + geom_vline(xintercept=0, linetype="dashed") + xlim(-0.5, 0.5) + scale_y_discrete(labels=c("AVG", "MOPAN(R)", "MOPAN(S)", "SWE", "NET", "DEN", "AUS", "UK")) + xlab("Coefficient estimate") + ylab(expression(Treatment:~Rating['i,a,t-1'])) + ggtitle("Baseline analysis") + theme(plot.title = element_text(hjust = 0.5)) + scale_color_manual(values=c("darkgray")) + geom_point(aes(x=rating_lgs$estimate, y=rating_lgs$term), colour="black")

## Eq. 4
## Fix conditioning variables at pre-treatment level
for(i in 1:length(assess)){
names <- paste(assess[i])
fix <- get(assess[i]) %>% group_by(IN) %>% mutate(comp_fix = ifelse(get(perf_vars[i]) != 0 & is.na(lag(get(perf_vars[i]), order_by=IN)), lag(competition, order_by=IN), NA), all_fix = ifelse(get(perf_vars[i]) != 0 & is.na(lag(get(perf_vars[i]), order_by=IN)), lag(alliances, order_by=IN), NA)) %>% fill(comp_fix, all_fix, .direction = "up") %>% fill(comp_fix, all_fix, .direction = "down")
assign(names, fix)
}

## Summary statistics
## TABLE A4 (Appendix)
data_stats <- as.data.frame(data[c("ln_cont_ia", "comp_fix", "all_fix", "rating_lg1")])
uk_stats <- as.data.frame(uk[c("rating_lg1", "comp_fix", "all_fix")])
aus_stats <- as.data.frame(aus[c("rating_lg1", "comp_fix", "all_fix")])
den_stats <- as.data.frame(den[c("rating_lg1", "comp_fix", "all_fix")])
net_stats <- as.data.frame(net[c("rating_lg1", "comp_fix", "all_fix")])
swe_stats <- as.data.frame(swe[c("rating_lg1", "comp_fix", "all_fix")])
mos_stats <- as.data.frame(mos[c("rating_lg1", "comp_fix", "all_fix")])
mor_stats <- as.data.frame(mor[c("rating_lg1", "comp_fix", "all_fix")])
stargazer(uk_stats, aus_stats, den_stats, net_stats, swe_stats, mos_stats, mor_stats, data_stats, type="html")

## Plotting moderator trends
## FIGURE A1 (Appendix)
pdf("Figure A1.pdf", width=15, height=8)
par(mfrow=c(1,2))
## Competition
comp_fig <- table(data$competition, data$year)
comp_fig <- comp_fig[,8:14]
barplot(comp_fig, xlab="Year", ylab="Frequency", main=expression(italic(Competition['i,a,t'])), beside=T, col=brewer.pal(5, "Blues"))
legend("topright", legend=rownames(comp_fig), fill=brewer.pal(5, "Blues"), horiz=T, xpd=T, inset=c(0,-0.025))
## Alliances
data$all_cats <- ifelse(data$alliances>=0 & data$alliances <= .25, 1, ifelse(data$alliances>.25 & data$alliances<=.5, 2, ifelse(data$alliances>.5 & data$alliances<=.75, 3, ifelse(data$alliances>.75, 4, 0))))
all_fig <- table(data$all_cats, data$year)
barplot(all_fig, xlab="Year", main=expression(italic(Alliances['i,a,t'])), beside=T, col=brewer.pal(6, "Reds"))
legend("top", legend=c("0-0.25", "0.25-0.5", "0.5-0.75", "0.75-1"), fill=brewer.pal(4, "Reds"), horiz=T, xpd=T, inset=F)
dev.off()

## Upper panel of TABLE 3
base4.lm <- base4.coef <- base4.se <- base4.p <- list()
for(i in 1:length(assess)){
	base4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess[i]))
	base4.coef[[i]] <- coeftest(base4.lm[[i]], cluster.vcov(base4.lm[[i]], get(assess[i])$IN))
	base4.se[[i]] <- base4.coef[[i]][,2]
	base4.p[[i]] <- base4.coef[[i]][,4]
}

stargazer(base4.lm, se=base4.se, p=base4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Marginal effects plots
## Factor
ip_factor <- function(model, effect, moderator1, interaction1, moderator2, interaction2, sample, varcov="default", conf=.95, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal effect", factor_labels=c(0,1,2,3,4), col1, col2){
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = cluster.vcov(model, sample$IN)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_4 = model$coefficients[[interaction1]]
  beta_5 = model$coefficients[[interaction2]]
 
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1,2,3,4)
  x_3 <- mean(mod_frame[[moderator2]], na.rm=T)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_4*x_2 + beta_5*x_3
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction1, interaction1] + (x_3^2)*covMat[interaction2, interaction2] + 2*x_2*covMat[effect, interaction1] + 2*x_3*covMat[effect, interaction2] + 2*x_2*x_3*covMat[interaction1, interaction2]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlab=xlabel, ylab=ylabel, xaxt="n", main=title, xlim=c(-0.5, 4.5))
  
  # Plot points of estimated effects
  points(x=x_2, y=delta_1, pch=16, col = col1)
  
  # Plot lines of confidence intervals
  lines(x=c(x_2[0], x_2[0]), y=c(upper_bound[0], lower_bound[0]), lty=1, col = col2)
  points(x=c(x_2[0], x_2[0]), y=c(upper_bound[0], lower_bound[0]), pch=c(25,24), bg= col2, col = col2)
  lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1, col = col2)
  points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg=col2, col = col2)
  lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1, col = col2)
  points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg=col2, col = col2)
    lines(x=c(x_2[3], x_2[3]), y=c(upper_bound[3], lower_bound[3]), lty=1, col = col2)
  points(x=c(x_2[3], x_2[3]), y=c(upper_bound[3], lower_bound[3]), pch=c(25,24), bg=col2, col = col2)
   lines(x=c(x_2[4], x_2[4]), y=c(upper_bound[4], lower_bound[4]), lty=1, col = col2)
  points(x=c(x_2[4], x_2[4]), y=c(upper_bound[4], lower_bound[4]), pch=c(25,24), bg=col2, col = col2)
     lines(x=c(x_2[5], x_2[5]), y=c(upper_bound[5], lower_bound[5]), lty=1, col = col2)
  points(x=c(x_2[5], x_2[5]), y=c(upper_bound[5], lower_bound[5]), pch=c(25,24), bg=col2, col = col2)
  
  # Label the axis
  axis(side=1, at=c(0,1,2,3,4), labels=factor_labels)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3, col = "gray50")
}

## Continuous
ip_continuous <- function(model, effect, moderator1, moderator2, interaction1, interaction2, sample, rugplot=T, varcov="default", minimum="min", maximum="max", incr="default", num_points = 50, conf=.95, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal effect", col1, col2){
    
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = cluster.vcov(model, sample$IN)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_4 = model$coefficients[[interaction1]]
  beta_5 = model$coefficients[[interaction2]]

  # Set range of the moderator variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[moderator1]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[moderator1]])
  }else{
    max_val = maximum
  }
  
  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum moderator value greater than maximum value.")
  }
  
  # Determine intervals between values of the moderator
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- seq(from=min_val, to=max_val, by=increment)
  x_3 <- mean(mod_frame[[moderator2]], na.rm=T)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_4*x_2 + beta_5*x_3
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction1, interaction1] + (x_3^2)*covMat[interaction2, interaction2] + 2*x_2*covMat[effect, interaction1] + 2*x_3*covMat[effect, interaction2] + 2*x_2*x_3*covMat[interaction1, interaction2]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title)
  
  # Plot estimated effects
  lines(y=delta_1, x=x_2, col = col1)
  lines(y=upper_bound, x=x_2, lty=2, col = col2)
  lines(y=lower_bound, x=x_2, lty=2, col = col2)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3, col = "gray50")
  
  ## Add rug plot
    if (rugplot){
    rug(mod_frame[[moderator1]])
  }
}

## MEs
base4.coef[[8]]
mean(data$all_fix)
mean(data$comp_fix)
quantile(data$all_fix, na.rm = T)
cor.test(data$comp_fix, data$all_fix)

## Unassessed control group
add <- read.csv("addins_data.csv", header=T)

## Log contributions
add$ln_cont_ia <- log(add$cont_ia+1)

## Exclude INs created after 2000
add_p2000 <- add[!(add$IN=="DCAF"|add$IN=="IDEA"|add$IN=="ITU"|add$IN=="WMO"),]
data_p2000 <- data[!(data$IN=="AF"|data$IN=="CERF"|data$IN=="CIFS"|data$IN=="EFW"|data$IN=="GCDT"|data$IN=="GPE"|data$IN=="GFATM"|data$IN=="GFDRR"|data$IN=="LDCF"|data$IN=="PBF"|data$IN=="PIDG"|data$IN=="UNITAID"|data$IN=="COMSEC"),]

## Comparison of funding trends
## Average values
add_mean <- aggregate(add_p2000[, "ln_cont_ia"], list(year=add_p2000 $year), mean, na.rm = T)
colnames(add_mean) <- c("year", "add_mean")
data_mean <- aggregate(data_p2000[, "ln_cont_ia"], list(year= data_p2000$year), mean, na.rm = T)
colnames(data_mean) <- c("year", "data_mean")
comparison <- merge(add_mean, data_mean, by="year")

## FIGURE 5
pdf("Figure 5.pdf")
legend_text <- c("Unassessed institutions", "Assessed Institutions", "First assessment")
plot_comp <- plot(comparison$year, comparison$add_mean, pch=17, xlab="Year", ylab="Log contributions (inflation-adjusted USD millions)", ylim=c(0, 10), main = NA, col ="black")
lines(comparison$year, comparison$add_mean, col = "black")
points(comparison$year, comparison$data_mean, pch=15, col = "darkgray")
lines(comparison$year, comparison$data_mean, lwd=0.6, col="darkgray")
abline(v=2008, lty=2, lwd=0.9, col="gray50")
legend("topleft", legend=legend_text, pch = c(17, 15, NA), lty=c(1, 1, 2), col=c("black", "darkgray", "gray50"))
dev.off()

## Add treatment and moderators
add$rating <- 0
add_uk <- add_aus <- add_den <- add_net <- add_swe <- add_mos <- add_mor <- add_data <- add
add_names <- c("add_uk", "add_aus", "add_den", "add_net", "add_swe", "add_mos", "add_mor", "add_data")
for(i in 1:length(add_names)){
	add$comp_fix <- mean(get(assess)$comp_fix, na.rm=T)
	add$all_fix <- mean(get(assess)$all_fix, na.rm=T)
	add <- add %>% group_by(IN) %>% mutate(rating_lg1 = lag(rating, order_by=IN))
	assign(paste(add_names[i]), add)
}

## Merge with existing data
for(i in 1:length(assess)){
	merged <- merge(add, get(assess[i]), all = T, by=c("IN", "year", "cont_ia", "ln_cont_ia", "rating_lg1", "comp_fix", "all_fix", "rating"))
	assign(paste(assess[i], "add", sep="_"), merged)
}

adds <- c("uk_add", "aus_add", "den_add", "net_add", "swe_add", "mos_add", "mor_add", "data_add")

## Leads and lags analysis
## Eq. 3
## TABLE A8 (Appendix)
lls3.lm <- lls3.coef <- lls3.se <- lls3.p <- list()
for(i in 1:length(assess)){
	lls3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating+rating_lg1+rating_lg2+rating_lg3+rating_ld1+rating_ld2+rating_ld3+factor(IN)+factor(year)")), data=get(adds[i]))
	lls3.coef[[i]] <- coeftest(lls3.lm[[i]], cluster.vcov(lls3.lm[[i]], get(adds[i])$IN))
	lls3.se[[i]] <- lls3.coef[[i]][,2]
	lls3.p[[i]] <- lls3.coef[[i]][,4]
}

stargazer(lls3.lm, se = lls3.se, p=lls3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 3
ucg3.lm <- ucg3.coef <- ucg3.se <- ucg3.p <- list()
for(i in 1:length(adds)){
	ucg3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)")), data=get(adds[i]))
	ucg3.coef[[i]] <- tidy(coeftest(ucg3.lm[[i]], cluster.vcov(ucg3.lm[[i]], get(adds[i])$IN)))
	ucg3.se[[i]] <- ucg3.coef[[i]][,2]
	ucg3.p[[i]] <- ucg3.coef[[i]][,4]
}

## Coefficient plot
l <- lapply(ucg3.coef, function(x) {i <- which(x$term == "rating_lg1")
	if (length(i) > 0) x[i, ]
	else NA })
rating_lgs <- rbind(l[[1]], l[[2]], l[[3]], l[[4]], l[[5]], l[[6]], l[[7]], l[[8]])
rating_lgs$term <- assess

coefplot2 <- dwplot(rating_lgs) + theme(panel.background = element_rect(fill = NA, colour = "black")) + geom_vline(xintercept=0, linetype="dashed") + xlim(-0.5, 0.5) + scale_y_discrete(labels=c("AVG", "MOPAN(R)", "MOPAN(S)", "SWE", "NET", "DEN", "AUS", "UK")) + xlab("Coefficient estimate") + ylab(expression(Treatment:~Rating['i,a,t-1'])) + ggtitle("Unassessed control group") + theme(plot.title = element_text(hjust = 0.5)) + scale_color_manual(values=c("darkgray")) + geom_point(aes(x=rating_lgs$estimate, y=rating_lgs$term), colour="black")

## Eq. 4
## Lower panel of TABLE 3
ucg4.lm <- ucg4.coef <- ucg4.se <- ucg4.p <- list()
for(i in 1:length(adds)){
	ucg4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(adds[i]))
	ucg4.coef[[i]] <- coeftest(ucg4.lm[[i]], cluster.vcov(ucg4.lm[[i]], get(adds[i])$IN))
	ucg4.se[[i]] <- ucg4.coef[[i]][,2]
	ucg4.p[[i]] <- ucg4.coef[[i]][,4]
}

stargazer(ucg4.lm, se=ucg4.se, p=ucg4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## ME plots
## FIGURE 4
pdf("Figure 4.pdf", width=10, height=10)
par(mfrow=c(2,2))
ip_factor(base4.lm[[8]], effect="rating_lg1", moderator1="comp_fix", interaction1="rating_lg1:comp_fix", moderator2="all_fix", interaction2="rating_lg1:all_fix", sample=data, factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~Rating['i,avg,t-1']), title = "Assessed Institutions Only", col1 = "black", col2="darkgray")
ip_continuous(base4.lm[[8]], effect="rating_lg1", moderator1="all_fix", moderator2="comp_fix", interaction1="rating_lg1:all_fix", interaction2="rating_lg1:comp_fix", sample=data, xlabel="Robustness of operational alliances", ylabel=NA, title = "Assessed Institutions Only", col1 = "black", col2="darkgray")
ip_factor(ucg4.lm[[8]], effect="rating_lg1", moderator1="comp_fix", interaction1="rating_lg1:comp_fix", moderator2="all_fix", interaction2="rating_lg1:all_fix", sample=data_add, factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~Rating['i,avg,t-1']), title = "Including Control Group", col1 = "black", col2="darkgray")
ip_continuous(ucg4.lm[[8]], effect="rating_lg1", moderator1="all_fix", moderator2="comp_fix", interaction1="rating_lg1:all_fix", interaction2="rating_lg1:comp_fix", sample= data_add, xlabel="Robustness of operational alliances", ylabel=NA, title = "Including Control Group", col1 = "black", col2="darkgray")
dev.off()

## RDD
## UK
## Running and treatment variables
uk$uk_sum <- ifelse(is.na(uk$uk_sum)==F, uk$uk_sum, 0)
uk$uk_sum2 <- uk$uk_sum^2
uk$uk_thresh <- ifelse(uk$uk_vfm==3 & !is.na(uk$uk_vfm), 1, ifelse(uk$uk_vfm==2& !is.na(uk$uk_vfm), -1, 0))

## Lags
uk <- uk %>% group_by(IN) %>% mutate(uk_sum_lg = lag(uk_sum, order_by=IN), uk_sum2_lg = lag(uk_sum2, order_by=IN), uk_thresh_lg = lag(uk_thresh, order_by=IN))

## Smallest bandwidth: uk_sum = [5, 6]
ukb56 <- uk %>% group_by(IN) %>% filter(sum(uk_sum)==5*6|sum(uk_sum)==6*6)
## Summary = 2 or 3: uk_sum = [5, 7]
ukb57 <- uk %>% group_by(IN) %>% filter(sum(uk_sum)>=5*6, sum(uk_sum)<=7*6)

## Eq. 5
ukbs <- c("ukb56", "ukb57")

rdd5_uk.lm <- rdd5_uk.coef <- rdd5_uk.se <- rdd5_uk.p <- list()
for(i in 1:length(ukbs)){
	rdd5_uk.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~uk_thresh_lg+uk_sum_lg+uk_sum2_lg+factor(IN)+factor(year)")), data=get(ukbs[i]))
	rdd5_uk.coef[[i]] <- tidy(coeftest(rdd5_uk.lm[[i]], cluster.vcov(rdd5_uk.lm[[i]], get(ukbs[i])$IN)))
	rdd5_uk.se[[i]] <- rdd5_uk.coef[[i]][,2]
	rdd5_uk.p[[i]] <- rdd5_uk.coef[[i]][,4]
}

## Eq. 6
rdd6_uk.lm <- rdd6_uk.coef <- rdd6_uk.se <- rdd6_uk.p <- list()
for(i in 1:length(ukbs)){
	rdd6_uk.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~uk_thresh_lg+ uk_thresh_lg:comp_fix+uk_thresh_lg:all_fix+uk_sum_lg+uk_sum2_lg+factor(IN)+factor(year)", sep="")), data=get(ukbs[i]))
	rdd6_uk.coef[[i]] <- coeftest(rdd6_uk.lm[[i]], cluster.vcov(rdd6_uk.lm[[i]], get(ukbs[i])$IN))
	rdd6_uk.se[[i]] <- rdd6_uk.coef[[i]][,2]
	rdd6_uk.p[[i]] <- rdd6_uk.coef[[i]][,4]
}

## Aus
## Running and treatment variables
aus$aus_int <- ifelse(is.na(aus$aus_int)==F, aus$aus_int, 0)
aus$aus_int2 <- aus$aus_int^2
aus$aus_thresh <- ifelse(aus$aus_vfm==3 & !is.na(aus$aus_vfm), 1, ifelse(aus$aus_vfm==2 & !is.na(aus$aus_vfm), -1, 0))

## Lags
aus <- aus %>% group_by(IN) %>% mutate(aus_int_lg = lag(aus_int, order_by=IN), aus_int2_lg = lag(aus_int2, order_by=IN), aus_thresh_lg = lag(aus_thresh, order_by=IN))

## Smallest band: aus_int = 3-4
ausb34 <- aus %>% group_by(IN) %>% filter(sum(aus_int)==3*5|sum(aus_int)==4*5)
## Summary = 2 or 3: aus_int = [2, 5]
ausb25 <- aus %>% group_by(IN) %>% filter(sum(aus_int)==2*5|sum(aus_int)==3*5|sum(aus_int)==4*5|sum(aus_int)==5*5)

## Eq. 5
ausbs <- c("ausb34", "ausb25")
rdd5_aus.lm <- rdd5_aus.coef <- rdd5_aus.se <- rdd5_aus.p <- list()
for(i in 1:length(ausbs)){
	rdd5_aus.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~aus_thresh_lg+aus_int_lg+aus_int2_lg+factor(IN)+factor(year)")), data=get(ausbs[i]))
	rdd5_aus.coef[[i]] <- tidy(coeftest(rdd5_aus.lm[[i]], cluster.vcov(rdd5_aus.lm[[i]], get(ausbs[i])$IN)))
	rdd5_aus.se[[i]] <- rdd5_aus.coef[[i]][,2]
	rdd5_aus.p[[i]] <- rdd5_aus.coef[[i]][,4]
}

## Eq. 6
rdd6_aus.lm <- rdd6_aus.coef <- rdd6_aus.se <- rdd6_aus.p <- list()
for(i in 1:length(ausbs)){
	rdd6_aus.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~aus_thresh_lg+aus_thresh_lg:comp_fix+aus_thresh_lg:all_fix+aus_int_lg+aus_int2_lg+factor(IN)+factor(year)", sep="")), data=get(ausbs[i]))
	rdd6_aus.coef[[i]] <- coeftest(rdd6_aus.lm[[i]], cluster.vcov(rdd6_aus.lm[[i]], get(ausbs[i])$IN))
	rdd6_aus.se[[i]] <- rdd6_aus.coef[[i]][,2]
	rdd6_aus.p[[i]] <- rdd6_aus.coef[[i]][,4]
}

## TABLE 4
stargazer(rdd6_uk.lm, se=rdd6_uk.se, p=rdd6_uk.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(rdd6_aus.lm, se=rdd6_aus.se, p=rdd6_aus.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Coefplot (Eq. 5)
l <- lapply(rdd5_uk.coef, function(x) {i <- which(x$term == "uk_thresh_lg")
                         if (length(i) > 0) x[i, ]
                         else NA })
rdduk <- rbind(l[[1]], l[[2]])
rdduk$term <- ukbs

l <- lapply(rdd5_aus.coef, function(x) {i <- which(x$term == "aus_thresh_lg")
                         if (length(i) > 0) x[i, ]
                         else NA })
rddaus <- rbind(l[[1]], l[[2]])
rddaus$term <- ausbs

rdd_rbind <- rbind(rdduk, rddaus)
coefplot3 <- dwplot(rdd_rbind) + theme(panel.background = element_rect(fill = NA, colour = "black")) + geom_vline(xintercept=0, linetype="dashed") + xlim(-1, 1) + scale_y_discrete(labels=c("AUS (band 2)", "AUS (band 1)", "UK (band 2)", "UK (band 1)")) + ylab(expression(Treatment:~Above/Below['i,a,t-1'])) + xlab("Coefficient estimate") + ggtitle("RDD") + theme(plot.title = element_text(hjust = 0.5)) + scale_color_manual(values=c("darkgray")) + geom_point(aes(x= rdd_rbind$estimate, y= rdd_rbind$term), colour="black")

## FIGURE 3
pdf("Figure 3.pdf", width=12, height=6)
grid.arrange(coefplot1, coefplot2, coefplot3, ncol=3)
dev.off()

## Marginal effects plots
## Band 1
## FIGURE A2
pdf("Figure A2.pdf", width=10, height=10)
par(mfrow=c(2,2))
ip_factor(rdd6_uk.lm[[2]], effect="uk_thresh_lg", moderator1="comp_fix", moderator2="all_fix", interaction1="uk_thresh_lg:comp_fix", interaction2="uk_thresh_lg:all_fix", sample=get(ukbs[2]), factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~Above/Below['i,uk,t-1']), col1="orange3", col2="orange", title="Band 2 (Running = [5, 7])")
ip_continuous(rdd6_uk.lm[[2]], effect="uk_thresh_lg", moderator1="all_fix", moderator2="comp_fix", interaction1="uk_thresh_lg:all_fix", interaction2="uk_thresh_lg:comp_fix", sample=get(ukbs[2]), xlabel="Robustness of operational alliances", ylabel=NA, col1="orange3", col2="orange", title ="Band 2 (Running = [5, 7])")
ip_factor(rdd6_aus.lm[[2]], effect="aus_thresh_lg", moderator1="comp_fix", moderator2="all_fix", interaction1="aus_thresh_lg:comp_fix", interaction2="aus_thresh_lg:all_fix", sample=get(ausbs[2]), factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~rating~Rating['i,aus,t-1']), title="Band 2 (Running = [2, 5])", col1="purple4", col2="purple")
ip_continuous(rdd6_aus.lm[[2]], effect="aus_thresh_lg", moderator1="all_fix", moderator2="comp_fix", interaction1="aus_thresh_lg:all_fix", interaction2="aus_thresh_lg:comp_fix", sample=get(ausbs[2]), xlabel="Robustness of operational alliances", ylabel=NA, title="Band 2 (Running = [2, 5])", col1="purple4", col2="purple")
dev.off()

## Band 2
## FIGURE A3
pdf("Figure A3.pdf", width=10, height=10)
par(mfrow=c(2,2))
ip_factor(rdd6_uk.lm[[1]], effect="uk_thresh_lg", moderator1="comp_fix", moderator2="all_fix", interaction1="uk_thresh_lg:comp_fix", interaction2="uk_thresh_lg:all_fix", sample=get(ukbs[1]), factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~Above/Below['i,uk,t-1']), col1="orange3", col2="orange", title="Band 1 (Running = [5, 6])")
ip_continuous(rdd6_uk.lm[[1]], effect="uk_thresh_lg", moderator1="all_fix", moderator2="comp_fix", interaction1="uk_thresh_lg:all_fix", interaction2="uk_thresh_lg:comp_fix", sample=get(ukbs[1]), xlabel="Robustness of operational alliances", ylabel=NA, col1="orange3", col2="orange", title ="Band 1 (Running = [5, 6])")
ip_factor(rdd6_aus.lm[[1]], effect="aus_thresh_lg", moderator1="comp_fix", moderator2="all_fix", interaction1="aus_thresh_lg:comp_fix", interaction2="aus_thresh_lg:all_fix", sample=get(ausbs[1]), factor_labels=c("0", "1", "2", "3", "4"), xlabel="Degree of institutional competition", ylabel=expression(Estimated~marginal~effect~of~rating~Rating['i,aus,t-1']), title="Band 1 (Running = [3, 4])", col1="purple4", col2="purple")
ip_continuous(rdd6_aus.lm[[1]], effect="aus_thresh_lg", moderator1="all_fix", moderator2="comp_fix", interaction1="aus_thresh_lg:all_fix", interaction2="aus_thresh_lg:comp_fix", sample=get(ausbs[1]), xlabel="Robustness of operational alliances", ylabel=NA, title="Band 1 (Running = [3, 4])", col1="purple4", col2="purple")
dev.off()

## McCrary test
## UK
DCdensity(uk$uk_sum, cutpoint=5.5)
## Aus
DCdensity(aus$aus_int, cutpoint=3.5, plot=F)
dev.off()

## Additional robustness checks
## (1) Additional controls
for(i in 1:length(assess)){
	vars <- cbind.data.frame(get(assess[i])$IN, get(assess[i])$year, get(assess[i])$rating)
	names(vars)[1:3] <- c("IN", "year", paste(assess[i], "_rating", sep=""))
	assign(paste(assess[i], "_vars", sep=""), vars)
}

merged <- list(data_vars, uk_vars, aus_vars, den_vars, net_vars, swe_vars, mos_vars, mor_vars) %>% Reduce(function(dtf1,dtf2) left_join(dtf1,dtf2,all=T,by=c("IN","year")), .)
merged$variance <- apply(merged[4:10],1,var,na.rm=T)
merged <- merged[,c(1,2,11)]
for(i in 1:length(assess)){
	merge <- merge(get(assess[i]),merged,by=c("IN", "year"))
	assign(paste(assess[i]), merge)
}

## Lag controls
for(i in 1:length(assess)){
	controls <- get(assess[i]) %>% group_by(IN) %>% mutate(var_lg = lag(variance, order_by=IN), vol_lg = lag(voluntary, order_by=IN), indep_lg = lag(independent, order_by=IN), ideal_lg = lag(ideal, order_by=IN))
	assign(paste(assess[i]), controls)
}

## Eq. 3
## TABLE A10 (Appendix)
cont3.lm <- cont3.coef <- cont3.se <- cont3.p <- list()
for(i in 1:length(assess)){
	cont3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+var_lg+vol_lg+indep_lg+ideal_lg+factor(IN)+factor(year)")), data=get(assess[i]))
	cont3.coef[[i]] <- coeftest(cont3.lm[[i]], cluster.vcov(cont3.lm[[i]], get(assess[i])$IN))
	cont3.se[[i]] <- cont3.coef[[i]][,2]
	cont3.p[[i]] <- cont3.coef[[i]][,4]
}

stargazer(cont3.lm, se=cont3.se, p=cont3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A11 (Appendix)
cont4.lm <- cont4.coef <- cont4.se <- cont4.p <- list()
for(i in 1:length(assess)){
	cont4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:comp_fix+rating_lg1:all_fix+var_lg+vol_lg+indep_lg+ideal_lg+factor(IN)+factor(year)", sep="")), data=get(assess[i]))
	cont4.coef[[i]] <- coeftest(cont4.lm[[i]], cluster.vcov(cont4.lm[[i]], get(assess[i])$IN))
	cont4.se[[i]] <- cont4.coef[[i]][,2]
	cont4.p[[i]] <- cont4.coef[[i]][,4]
}

stargazer(cont4.lm, se=cont4.se, p=cont4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## (2) Aggregation to two periods 
for(i in 1:length(assess)){
	post <- get(assess[i])[which(!get(assess[i])$rating==0),]
	post$period=1
	pre <- get(assess[i])[which(get(assess[i])$rating==0),]
	pre$period=2
	post <- aggregate(post[, c("ln_cont_ia", "comp_fix", "all_fix", "rating_lg1", "period")], list(IN=post$IN), mean, na.rm = T)
	pre <- aggregate(pre[, c("ln_cont_ia", "comp_fix", "all_fix", "rating_lg1", "period")], list(IN=pre$IN), mean, na.rm = T)
	agg <- rbind(post, pre)
	assign(paste(assess[i], "agg", sep="_"), agg)
}

aggs <- c("uk_agg", "aus_agg", "den_agg", "net_agg", "swe_agg", "mos_agg", "mor_agg", "data_agg")

## Eq. 3
## TABLE A12 (Appendix)
per3.lm <- per3.coef <- per3.se <- per3.p <- list()
for(i in 1:length(aggs)){
	per3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(period)")), data=get(aggs[i]))
	per3.coef[[i]] <- coeftest(per3.lm[[i]], cluster.vcov(per3.lm[[i]], get(aggs[i])$IN))
	per3.se[[i]] <- per3.coef[[i]][,2]
	per3.p[[i]] <- per3.coef[[i]][,4]
}

stargazer(per3.lm, se=per3.se, p=per3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"), omit.stat=c("ser","f"))

## Eq. 4
## TABLE A13 (Appendix)
per4.lm <- per4.coef <- per4.se <- per4.p <- list()
for(i in 1:length(aggs)){
	per4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(period)", sep="")), data=get(aggs[i]))
	per4.coef[[i]] <- coeftest(per4.lm[[i]], cluster.vcov(per4.lm[[i]], get(aggs[i])$IN))
	per4.se[[i]] <- per4.coef[[i]][,2]
	per4.p[[i]] <- per4.coef[[i]][,4]
}

stargazer(per4.lm, se=per4.se, p=per4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## (3) Separating moderators and interactions
## comp_fix only
## TABLE A14 (Appendix)
sep_comp4.lm <- sep_comp4.coef <- sep_comp4.se <- sep_comp4.p <- list()
for(i in 1:length(assess)){
	sep_comp4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:comp_fix+factor(IN)+factor(year)", sep_comp="")), data=get(assess[i]))
	sep_comp4.coef[[i]] <- coeftest(sep_comp4.lm[[i]], cluster.vcov(sep_comp4.lm[[i]], get(assess[i])$IN))
	sep_comp4.se[[i]] <- sep_comp4.coef[[i]][,2]
	sep_comp4.p[[i]] <- sep_comp4.coef[[i]][,4]
}

stargazer(sep_comp4.lm, se=sep_comp4.se, p=sep_comp4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## all_fix only
## TABLE A15 (Appendix)
sep_all4.lm <- sep_all4.coef <- sep_all4.se <- sep_all4.p <- list()
for(i in 1:length(assess)){
	sep_all4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:all_fix+factor(IN)+factor(year)", sep_all="")), data=get(assess[i]))
	sep_all4.coef[[i]] <- coeftest(sep_all4.lm[[i]], cluster.vcov(sep_all4.lm[[i]], get(assess[i])$IN))
	sep_all4.se[[i]] <- sep_all4.coef[[i]][,2]
	sep_all4.p[[i]] <- sep_all4.coef[[i]][,4]
}
	
stargazer(sep_all4.lm, se=sep_all4.se, p=sep_all4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## (4) Triple interaction term
## TABLE A16 (Appendix)
triple4.lm <- triple4.coef <- triple4.se <- triple4.p <- list()
for(i in 1:length(assess)){
	triple4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:comp_fix + rating_lg1:all_fix+comp_fix:all_fix+rating_lg1:comp_fix:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess[i]))
	triple4.coef[[i]] <- coeftest(triple4.lm[[i]], cluster.vcov(triple4.lm[[i]], get(assess[i])$IN))
	triple4.se[[i]] <- triple4.coef[[i]][,2]
	triple4.p[[i]] <- triple4.coef[[i]][,4]
}

stargazer(triple4.lm, se= triple4.se, p= triple4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## (5) Alternative interaction estimation strategies
## (a) Binning estimates
data_df <- as.data.frame(data)

## FIGURE A4 (Appendix)
## Competition
pdf("Figure A4.pdf", width=14, height=7)
comp_bin <- inter.binning(Y = "ln_cont_ia", D = "rating_lg1", X = "comp_fix", data = data_df, nbins = 5, FE = c("IN", "year"), cl = "IN", na.rm=T, Ylabel="Log Contributions (t)", Xlabel="Competition (g-1)", Dlabel="Avg. Rating (t-1)")
## Alliances
all_bin <- inter.binning(Y = "ln_cont_ia", D = "rating_lg1", X = "all_fix", data = data_df, nbins = 5, FE = c("IN", "year"), cl = "IN", na.rm=T, Ylabel="Log Contributions (t)", Xlabel="Alliances (g-1)", Dlabel="Avg. Rating (t-1)")
grid.arrange(comp_bin$graph, all_bin$graph, ncol=2)
dev.off()

## (b) Kernel estimates
## FIGURE A5 (Appendix)
## Competition
pdf("Figure A5.pdf", width=14, height=7)
comp_ker <- inter.kernel(Y = "ln_cont_ia", D = "rating_lg1", X = "comp_fix", data = data_df, FE = c("IN", "year"), cl = "IN", na.rm=T, Ylabel="Log Contributions (t)", Xlabel="Competition (g-1)", Dlabel="Avg. Rating (t-1)", weights=NULL)
## Alliances
all_ker <- inter.kernel(Y = "ln_cont_ia", D = "rating_lg1", X = "all_fix", data = data_df, FE = c("IN", "year"), cl = "IN", na.rm=T, Ylabel="Log Contributions (t)", Xlabel="Alliances (g-1)", Dlabel="Avg. Rating (t-1)", weights=NULL)
grid.arrange(comp_ker$graph, all_ker$graph, ncol=2)
dev.off()

## (6) Alternative measure of competition
## Fix at pre-treatment level
for(i in 1:length(assess)){
fix <- get(assess[i]) %>% group_by(IN) %>% mutate(nstand_fix = ifelse(get(perf_vars[i]) != 0 & is.na(lag(get(perf_vars[i]), order_by=IN)), lag(1-standards, order_by=IN), NA)) %>% fill(nstand_fix, .direction = "up") %>% fill(nstand_fix, .direction = "down")
assign(paste(assess[i]), fix)
}

## TABLE A17 (Appendix)
nstand4.lm <- nstand4.coef <- nstand4.se <- nstand4.p <- list()
for(i in 1:length(assess)){
	nstand4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+ rating_lg1:nstand_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess[i]))
	nstand4.coef[[i]] <- coeftest(nstand4.lm[[i]], cluster.vcov(nstand4.lm[[i]], get(assess[i])$IN))
	nstand4.se[[i]] <- nstand4.coef[[i]][,2]
	nstand4.p[[i]] <- nstand4.coef[[i]][,4]
}

stargazer(nstand4.lm, se=nstand4.se, p=nstand4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Interactive fixed effects
## Development
## Eq. 3
## TABLE A18 (Appendix)
intfe_dev3.lm <- intfe_dev3.coef <- intfe_dev3.se <- intfe_dev3.p <- list()
for(i in 1:length(assess)){
	intfe_dev3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)*development")), data=get(assess[i]))
	intfe_dev3.coef[[i]] <- coeftest(intfe_dev3.lm[[i]], cluster.vcov(intfe_dev3.lm[[i]], get(assess[i])$IN))
	intfe_dev3.se[[i]] <- intfe_dev3.coef[[i]][,2]
	intfe_dev3.p[[i]] <- intfe_dev3.coef[[i]][,4]
}

stargazer(intfe_dev3.lm, se= intfe_dev3.se, p= intfe_dev3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A19 (Appendix)
intfe_dev4.lm <- intfe_dev4.coef <- intfe_dev4.se <- intfe_dev4.p <- list()
for(i in 1:length(assess)){
	intfe_dev4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)*development", sep="")), data=get(assess[i]))
	intfe_dev4.coef[[i]] <- coeftest(intfe_dev4.lm[[i]], cluster.vcov(intfe_dev4.lm[[i]], get(assess[i])$IN))
	intfe_dev4.se[[i]] <- intfe_dev4.coef[[i]][,2]
	intfe_dev4.p[[i]] <- intfe_dev4.coef[[i]][,4]
}

stargazer(intfe_dev4.lm, se= intfe_dev4.se, p= intfe_dev4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Education
## Eq. 3
## TABLE A20 (Appendix)
intfe_edu3.lm <- intfe_edu3.coef <- intfe_edu3.se <- intfe_edu3.p <- list()
for(i in 1:length(assess)){
	intfe_edu3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)*education")), data=get(assess[i]))
	intfe_edu3.coef[[i]] <- coeftest(intfe_edu3.lm[[i]], cluster.vcov(intfe_edu3.lm[[i]], get(assess[i])$IN))
	intfe_edu3.se[[i]] <- intfe_edu3.coef[[i]][,2]
	intfe_edu3.p[[i]] <- intfe_edu3.coef[[i]][,4]
}

stargazer(intfe_edu3.lm, se= intfe_edu3.se, p= intfe_edu3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A21 (Appendix)
intfe_edu4.lm <- intfe_edu4.coef <- intfe_edu4.se <- intfe_edu4.p <- list()
for(i in 1:length(assess)){
	intfe_edu4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)*education", sep="")), data=get(assess[i]))
	intfe_edu4.coef[[i]] <- coeftest(intfe_edu4.lm[[i]], cluster.vcov(intfe_edu4.lm[[i]], get(assess[i])$IN))
	intfe_edu4.se[[i]] <- intfe_edu4.coef[[i]][,2]
	intfe_edu4.p[[i]] <- intfe_edu4.coef[[i]][,4]
}

stargazer(intfe_edu4.lm, se= intfe_edu4.se, p= intfe_edu4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Environment
## Eq. 3
## TABLE A22 (Appendix)
intfe_env3.lm <- intfe_env3.coef <- intfe_env3.se <- intfe_env3.p <- list()
for(i in 1:length(assess)){
	intfe_env3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)*environment")), data=get(assess[i]))
	intfe_env3.coef[[i]] <- coeftest(intfe_env3.lm[[i]], cluster.vcov(intfe_env3.lm[[i]], get(assess[i])$IN))
	intfe_env3.se[[i]] <- intfe_env3.coef[[i]][,2]
	intfe_env3.p[[i]] <- intfe_env3.coef[[i]][,4]
}

stargazer(intfe_env3.lm, se= intfe_env3.se, p= intfe_env3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A23 (Appendix)
intfe_env4.lm <- intfe_env4.coef <- intfe_env4.se <- intfe_env4.p <- list()
for(i in 1:length(assess)){
	intfe_env4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)*environment", sep="")), data=get(assess[i]))
	intfe_env4.coef[[i]] <- coeftest(intfe_env4.lm[[i]], cluster.vcov(intfe_env4.lm[[i]], get(assess[i])$IN))
	intfe_env4.se[[i]] <- intfe_env4.coef[[i]][,2]
	intfe_env4.p[[i]] <- intfe_env4.coef[[i]][,4]
}

stargazer(intfe_env4.lm, se= intfe_env4.se, p= intfe_env4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Humanitarian
## Eq. 3
## TABLE A24 (Appendix)
intfe_hum3.lm <- intfe_hum3.coef <- intfe_hum3.se <- intfe_hum3.p <- list()
for(i in 1:length(assess)){
	intfe_hum3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)*humanitarian")), data=get(assess[i]))
	intfe_hum3.coef[[i]] <- coeftest(intfe_hum3.lm[[i]], cluster.vcov(intfe_hum3.lm[[i]], get(assess[i])$IN))
	intfe_hum3.se[[i]] <- intfe_hum3.coef[[i]][,2]
	intfe_hum3.p[[i]] <- intfe_hum3.coef[[i]][,4]
}

stargazer(intfe_hum3.lm, se= intfe_hum3.se, p= intfe_hum3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A25 (Appendix)
intfe_hum4.lm <- intfe_hum4.coef <- intfe_hum4.se <- intfe_hum4.p <- list()
for(i in 1:length(assess)){
	intfe_hum4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)*humanitarian", sep="")), data=get(assess[i]))
	intfe_hum4.coef[[i]] <- coeftest(intfe_hum4.lm[[i]], cluster.vcov(intfe_hum4.lm[[i]], get(assess[i])$IN))
	intfe_hum4.se[[i]] <- intfe_hum4.coef[[i]][,2]
	intfe_hum4.p[[i]] <- intfe_hum4.coef[[i]][,4]
}

stargazer(intfe_hum4.lm, se= intfe_hum4.se, p= intfe_hum4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Health
## Eq. 3
## TABLE A26 (Appendix)
intfe_hea3.lm <- intfe_hea3.coef <- intfe_hea3.se <- intfe_hea3.p <- list()
for(i in 1:length(assess)){
	intfe_hea3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)*health")), data=get(assess[i]))
	intfe_hea3.coef[[i]] <- coeftest(intfe_hea3.lm[[i]], cluster.vcov(intfe_hea3.lm[[i]], get(assess[i])$IN))
	intfe_hea3.se[[i]] <- intfe_hea3.coef[[i]][,2]
	intfe_hea3.p[[i]] <- intfe_hea3.coef[[i]][,4]
}

stargazer(intfe_hea3.lm, se= intfe_hea3.se, p= intfe_hea3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A27 (Appendix)
intfe_hea4.lm <- intfe_hea4.coef <- intfe_hea4.se <- intfe_hea4.p <- list()
for(i in 1:length(assess)){
	intfe_hea4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)*health", sep="")), data=get(assess[i]))
	intfe_hea4.coef[[i]] <- coeftest(intfe_hea4.lm[[i]], cluster.vcov(intfe_hea4.lm[[i]], get(assess[i])$IN))
	intfe_hea4.se[[i]] <- intfe_hea4.coef[[i]][,2]
	intfe_hea4.p[[i]] <- intfe_hea4.coef[[i]][,4]
}

stargazer(intfe_hea4.lm, se= intfe_hea4.se, p= intfe_hea4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN"), omit.labels=c("Institution FE"),omit.stat=c("ser","f"))

## UN vs. other institutions
## Lag and subset data
for(i in 1:length(assess)){
	assign(paste0(assess[i], "_un"), get(assess[i])[ get(assess[i])$un==1, ])
	assign(paste0(assess[i], "_nonun"), get(assess[i])[ get(assess[i])$un==0, ])
}

assess_un <- paste0(assess, "_un")
assess_nonun <- paste0(assess, "_nonun")

## Interaction term
## Eq. 3
## TABLE A28 (Appendix)
base_unint3.lm <- base_unint3.coef <- base_unint3.se <- base_unint3.p <- list()
for(i in 1:length(assess)){
	base_unint3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:un+factor(IN)+factor(year)")), data=get(assess[i]))
	base_unint3.coef[[i]] <- coeftest(base_unint3.lm[[i]], cluster.vcov(base_unint3.lm[[i]], get(assess[i])$IN))
	base_unint3.se[[i]] <- base_unint3.coef[[i]][,2]
	base_unint3.p[[i]] <- base_unint3.coef[[i]][,4]
}

stargazer(base_unint3.lm, se= base_unint3.se, p= base_unint3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 4
## TABLE A29 (Appendix)
base_unint4.lm <- base_unint4.coef <- base_unint4.se <- base_unint4.p <- list()
for(i in 1:length(assess)){
	base_unint4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+rating_lg1:un+rating_lg1:un:comp_fix+rating_lg1:un:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess[i]))
	base_unint4.coef[[i]] <- coeftest(base_unint4.lm[[i]], cluster.vcov(base_unint4.lm[[i]], get(assess[i])$IN))
	base_unint4.se[[i]] <- base_unint4.coef[[i]][,2]
	base_unint4.p[[i]] <- base_unint4.coef[[i]][,4]
}

stargazer(base_unint4.lm, se= base_unint4.se, p= base_unint4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Split sample
## Eq. 3: UN
## TABLE A30 (Appendix)
base_un3.lm <- base_un3.coef <- base_un3.se <- base_un3.p <- list()
for(i in 1:length(assess)){
	base_un3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)")), data=get(assess_un[i]))
	base_un3.coef[[i]] <- coeftest(base_un3.lm[[i]], cluster.vcov(base_un3.lm[[i]], get(assess_un[i])$IN))
	base_un3.se[[i]] <- base_un3.coef[[i]][,2]
	base_un3.p[[i]] <- base_un3.coef[[i]][,4]
}

stargazer(base_un3.lm, se= base_un3.se, p= base_un3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 4: UN
## TABLE A31 (Appendix)
base_un4.lm <- base_un4.coef <- base_un4.se <- base_un4.p <- list()
for(i in 1:length(assess)){
	base_un4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess_un[i]))
	base_un4.coef[[i]] <- coeftest(base_un4.lm[[i]], cluster.vcov(base_un4.lm[[i]], get(assess_un[i])$IN))
	base_un4.se[[i]] <- base_un4.coef[[i]][,2]
	base_un4.p[[i]] <- base_un4.coef[[i]][,4]
}

stargazer(base_un4.lm, se= base_un4.se, p= base_un4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 3: Non-UN
## TABLE A32 (Appendix)
base_nonun3.lm <- base_nonun3.coef <- base_nonun3.se <- base_nonun3.p <- list()
for(i in 1:length(assess)){
	base_nonun3.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+factor(IN)+factor(year)")), data=get(assess_nonun[i]))
	base_nonun3.coef[[i]] <- coeftest(base_nonun3.lm[[i]], cluster.vcov(base_nonun3.lm[[i]], get(assess_nonun[i])$IN))
	base_nonun3.se[[i]] <- base_nonun3.coef[[i]][,2]
	base_nonun3.p[[i]] <- base_nonun3.coef[[i]][,4]
}

stargazer(base_nonun3.lm, se= base_nonun3.se, p= base_nonun3.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Eq. 4: Non-UN
## TABLE A33 (Appendix)
base_nonun4.lm <- base_nonun4.coef <- base_nonun4.se <- base_nonun4.p <- list()
for(i in 1:length(assess)){
	base_nonun4.lm[[i]] <- lm(as.formula(paste("ln_cont_ia~rating_lg1+rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess_nonun[i]))
	base_nonun4.coef[[i]] <- coeftest(base_nonun4.lm[[i]], cluster.vcov(base_nonun4.lm[[i]], get(assess_nonun[i])$IN))
	base_nonun4.se[[i]] <- base_nonun4.coef[[i]][,2]
	base_nonun4.p[[i]] <- base_nonun4.coef[[i]][,4]
}

stargazer(base_nonun4.lm, se= base_nonun4.se, p= base_nonun4.p, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Disaggregation of ratings and contributions
## Sub-indicator analysis
## Create sub-indicator vectors
sub_uk <- c("uk_clim", "uk_res", "uk_cost", "uk_critint", "uk_critnat", "uk_fin", "uk_poor", "uk_frag", "uk_gend", "uk_part", "uk_strat", "uk_trans")
sub_aus <- c("aus_natint", "aus_multi", "aus_cost", "aus_res", "aus_part", "aus_strat", "aus_trans")
sub_den <- c("den_align", "den_comp", "den_cont", "den_fin", "den_actors", "den_innov", "den_reform", "den_part", "den_info")
sub_net <- c("net_cor", "net_gov", "net_fin", "net_hrm", "net_part", "net_eval", "net_res", "net_strat")
sub_swe <- c("swe_ext", "swe_int", "swe_relev")
sub_mos <- c("mos_AP", "mos_CPD", "mos_CTFR", "mos_DLL", "mos_EXR", "mos_FA", "mos_FTP", "mos_HP", "mos_PDR", "mos_RAD", "mos_UPI")
sub_mor <- c("mor_CPFR", "mor_DLL", "mor_EXR", "mor_FA", "mor_FTP", "mor_MHR", "mor_PPI", "mor_PDR", "mor_RAD", "mor_UPI")

## Define treatments 
subs <- list(sub_uk, sub_aus, sub_den, sub_net, sub_swe, sub_mos, sub_mor)
assess2 <- head(assess, -1)
res_rating <- names_rating <- names_rating_lg <- vector('list', length(subs))
for(i in seq_along(subs)){
	for(j in seq_along(subs[[i]])){
	names_rating[[i]][[j]] <- paste(subs[[i]][[j]], "_rating", sep="")
	res_rating[[i]][[j]] <- ifelse(!is.na(get(assess2[i])[[subs[[i]][[j]]]]),scale(get(assess2[i])[[subs[[i]][[j]]]]), 0)
	names(res_rating)[[i]] <- assess2[[i]]
	names(res_rating[[i]]) <- names_rating[[i]]
	unlist_rating <- as.data.frame(do.call(cbind, res_rating[[i]]))
	bind <- bind_cols(get(assess2[i]), unlist_rating)
	assign(paste(assess2[[i]]), bind)
	}
}

## Lags
names_rating_lg <- vector('list', length(subs))
for(i in 1:length(assess2)){
lags <- get(assess2[i]) %>% group_by(IN) %>% mutate_at(vars(names_rating[[i]]), list("lg"=lag)) %>% ungroup
	assign(paste(assess2[[i]]), lags)
	names_rating_lg[[i]] <- paste(names_rating[[i]], "_lg", sep="")
}

## Eq. 3
sub3.lm <- sub3.coef <- sub3.se <- sub3.p <- vector('list', length(subs))
for(i in seq_along(subs)){
	for(j in seq_along(subs[[i]])){
	sub3.lm[[i]][[j]] <- lm(as.formula(paste("ln_cont_ia~",names_rating_lg[[i]][[j]],"+factor(IN)+factor(year)")), data=get(assess[i]))
	sub3.coef[[i]][[j]] <- coeftest(sub3.lm[[i]][[j]], cluster.vcov(sub3.lm[[i]][[j]], get(assess[i])$IN))
	sub3.se[[i]][[j]] <- sub3.coef[[i]][[j]][,2]
	sub3.p[[i]][[j]] <- sub3.coef[[i]][[j]][,4]
		}
}

## Eq. 4
sub4.lm <- sub4.coef <- sub4.se <- sub4.p <- vector('list', length(subs))
for(i in seq_along(subs)){
	for(j in seq_along(subs[[i]])){
	sub4.lm[[i]][[j]] <- lm(as.formula(paste("ln_cont_ia~",names_rating_lg[[i]][[j]],"*comp_fix+",names_rating_lg[[i]][[j]],"*all_fix+factor(IN)+factor(year)")), data=get(assess2[i]))
	sub4.coef[[i]][[j]] <- coeftest(sub4.lm[[i]][[j]], cluster.vcov(sub4.lm[[i]][[j]], get(assess[i])$IN))
	sub4.se[[i]][[j]] <- sub4.coef[[i]][[j]][,2]
	sub4.p[[i]][[j]] <- sub4.coef[[i]][[j]][,4]
	}
}

## UK
## TABLE A35 (Appendix)
stargazer(sub3.lm[[1]], se=sub3.se[[1]], p=sub3.p[[1]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[1]], se=sub4.se[[1]], p=sub4.p[[1]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Aus
## TABLE A36 (Appendix)
stargazer(sub3.lm[[2]], se=sub3.se[[2]], p=sub3.p[[2]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[2]], se=sub4.se[[2]], p=sub4.p[[2]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Den
## TABLE A37 (Appendix)
stargazer(sub3.lm[[3]], se=sub3.se[[3]], p=sub3.p[[3]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[3]], se=sub4.se[[3]], p=sub4.p[[3]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Net
## TABLE A38 (Appendix)
stargazer(sub3.lm[[4]], se=sub3.se[[4]], p=sub3.p[[4]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[4]], se=sub4.se[[4]], p=sub4.p[[4]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Swe
## TABLE A39 (Appendix)
stargazer(sub3.lm[[5]], se=sub3.se[[5]], p=sub3.p[[5]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[5]], se=sub4.se[[5]], p=sub4.p[[5]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Mos
## TABLE A40 (Appendix)
stargazer(sub3.lm[[6]], se=sub3.se[[6]], p=sub3.p[[6]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[6]], se=sub4.se[[6]], p=sub4.p[[6]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Mor
## TABLE A41 (Appendix)
stargazer(sub3.lm[[7]], se=sub3.se[[7]], p=sub3.p[[7]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))
stargazer(sub4.lm[[7]], se=sub4.se[[7]], p=sub4.p[[7]], align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Analysis of individual donors' contributions
nat_cont <- read.csv("natcont_data.csv", header=T)
indiv <- c("aus", "bel", "can", "den", "fin", "fra", "ger", "ire", "ita", "jap", "kor", "lux", "net", "nor", "swe", "swi", "uk", "usa")
for(i in 1:length(indiv)){
	nat_cont[[paste("ln",indiv[i],"ia", sep="_")]] <-ifelse(nat_cont[[indiv[i]]]<0, 0, log(nat_cont[[indiv[i]]]+1))
}

## Merge
for(i in 1:length(assess)){
	merged <- merge(get(assess[i]), nat_cont, all=T, by=c("IN", "year"))
	assign(paste(assess[i]), merged)
}

## Eq. 3
indiv3.lm <- indiv3.coef <- indiv3.se <- indiv3.p <- vector('list', length(indiv))
for(i in 1:length(indiv)){
	for(j in 1:length(assess)){
	indiv3.lm[[i]][[j]] <- lm(as.formula(paste("ln_",indiv[i],"_ia","~rating_lg1+factor(IN)+factor(year)", sep="")), data=get(assess[j]))
	indiv3.coef[[i]][[j]] <- coeftest(indiv3.lm[[i]][[j]], cluster.vcov(indiv3.lm[[i]][[j]], get(assess[j])$IN))
	indiv3.se[[i]][[j]] <- indiv3.coef[[i]][[j]][,2]
	indiv3.p[[i]][[j]] <- indiv3.coef[[i]][[j]][,4]
	}
}

## Eq. 4
indiv4.lm <- indiv4.coef <- indiv4.se <- indiv4.p <- vector('list', length(indiv))
for(i in 1:length(indiv)){
	for(j in 1:length(assess)){
	indiv4.lm[[i]][[j]] <- lm(as.formula(paste("ln_",indiv[i],"_ia","~rating_lg1+ rating_lg1:comp_fix+rating_lg1:all_fix+factor(IN)+factor(year)", sep="")), data=get(assess[j]))
	indiv4.coef[[i]][[j]] <- coeftest(indiv4.lm[[i]][[j]], cluster.vcov(indiv4.lm[[i]][[j]], get(assess[j])$IN))
	indiv4.se[[i]][[j]] <- indiv4.coef[[i]][[j]][,2]
	indiv4.p[[i]][[j]] <- indiv4.coef[[i]][[j]][,4]
	}
}

## Summary
for(i in 1:length(indiv)){
	for(j in 1:length(assess)){
	rating3 <- indiv3.coef[[i]][[j]]["rating_lg1",c(1,2,4)]
	ratingcomp4 <- indiv4.coef[[i]][[j]]["rating_lg1:comp_fix",c(1,2,4)]
	ratingall4 <- indiv4.coef[[i]][[j]]["rating_lg1:all_fix",c(1,2,4)]
	combined <- cbind(rating3, ratingcomp4, ratingall4)
	assign(paste("sum",i,j, sep="."), combined)
	}
}

## Call results by assessment
for(i in 1:length(assess)){
	res <- do.call(rbind, lapply( paste("sum",1:18,i,sep=".") , get) )
	assign(paste(assess[i],"res",sep="_"),res)
}

## TABLE A42 (Appendix)
stargazer(uk_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A43 (Appendix)
stargazer(aus_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A44 (Appendix)
stargazer(den_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A45 (Appendix)
stargazer(net_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A46 (Appendix)
stargazer(swe_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A47 (Appendix)
stargazer(mos_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A48 (Appendix)
stargazer(mor_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## TABLE A49 (Appendix)
stargazer(data_res, align = TRUE, type = "html", digits = 3, no.space = T, omit=c("IN","year"), omit.labels=c("Institution FE", "Year FE"),omit.stat=c("ser","f"))

## Call results by donor
for(i in 1:length(indiv)){
	ind_res <- do.call(rbind, lapply( paste("sum", i, 1:8, sep=".") , get) )
	assign(paste(indiv[i],"ind_res",sep="_"), ind_res)
}